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BLOCK-DIAGONALIZATION OF THE LINEARIZED 
COUPLED-MODE SYSTEM 

MARINA CHUGUNOVA, DMITRY PELINOVSKY 

Abstract. We consider the Hamiltonian coupled-mode system that occur 
in nonlinear optics, photonics, and atomic physics. Spectral stability of gap 
solitons is determined by eigenvalues of the linearized coupled-mode system, 
which is equivalent to a four-by-four Dirac system with sign-indefinite metric. 
In the special class of symmetric nonlinear potentials, we construct a block- 
diagonal representation of the linearized equations, when the spectral problem 
reduces to two coupled two-by-two Dirac systems. The block-diagonalization 
is used in numerical computations of eigenvalues that determine stability of 
gap solitons. 



1. Introduction 

Various applications in nonlinear optics 1 , photonics band-gap engineering |2] 
and atomic physics j2| call for systematic studies of the coupled-mode system, which 
is expressed by two first-order semi-linear PDEs in one space and one time dimen- 
sions. In nonlinear optics, the coupled-mode system describes counter-propagating 
light waves, which interact with a linear grating in an optical waveguide 0]. In 
photonics, the coupled-mode system is derived for coupled resonant waves in stop 
bands of a low-contrast three-dimensional photonic crystal [5]. In atomic physics, 
the coupled-mode system describes matter- wave Bose-Einstein condensates trapped 
in an optical lattice W. Existence, stability and nonlinear dynamics of gap soli- 
tons, which are localized solutions of the coupled-mode system, are fundamental 
problems for interest in the aforementioned physical disciplines. 

In the context of spectral stability of gap solitons, it has been discovered that 
the linearized coupled-mode system is equivalent to a four-by-four Dirac system 
with sign-indefinite metric, where numerical computations of eigenvalues represent 
a difficult numerical task. The pioneer work in j3|8j showed that spurious unstable 
eigenvalues originate from the continuous spectrum in the Fourier basis decompo- 
sition and the Galerkin approximation. A delicate but time-consuming implemen- 
tation of the continuous Newton method was developed to identify true unstable 
eigenvalues from the spurious ones [71 |S] . Similar problems were discovered in the 
variational method [HI EI] and in the numerical finite-difference method ^] E] • 

While some conclusions on instability bifurcations of gap solitons in the coupled- 
mode equations can be drawn on the basis of perturbation theory and Evans 
function methods [131 114| , the numerical approximation of eigenvalues was an open 
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problem until recently. A new progress was made with the use of exterior algebra 
in the numerical computations of the Evans function [ISj , when the same results 
on instability bifurcations of gap solitons as in were recovered. Similar shooting 
method was also applied to gap solitons in a more general model of a nonlinear 
Schrodinger equation with a periodic potential 

Our work addresses the problem of numerical approximations of eigenvalues of 
the linearized coupled-mode system with a different objective. We will show that 
the linearized coupled-mode system with a symmetric potential function can be 
block-diagonalized into two coupled two-by-two Dirac systems. The two Dirac 
systems represent the linearized Hamiltonian of the coupled-mode equations and 
determine instability bifurcations and unstable eigenvalues of gap solitons. 

The purpose of block-diagonalization is twofold. First, the number of unstable 
eigenvalues and details of instability bifurcations can be investigated analytically 
from the number of non-zero isolated eigenvalues of the linearized Hamiltonian. 
This analysis will be reported in the forthcoming publication. Second, a numeri- 
cal algorithm can be developed to compute efficiently the entire spectrum of the 
linearized coupled-mode system. These numerical results are reported here for an 
example of symmetric quadric potential functions. 

The paper is organized as follows. Section 2 describes the model and its symme- 
tries. Section 3 gives construction and properties of gap solitons in the nonlinear 
coupled-mode system. Section 4 presents block-diagonalization of the linearized 
coupled-mode system. Section 5 contains numerical computations of the spectrum 
of the block-diagonalized system. Appendix A presents derivation of exact solu- 
tions for gap solitons in the coupled-mode system with symmetric quadric potential 
functions. 



where {u,v) £ C'^, x E M., t > 0, and W{u,u,v,v) is real-valued. We assume that 
the potential function satisfies the following three conditions: 

(1) W is invariant with respect to the gauge transformation: (u, v) i— > e'"(u, v), 
for all a G E 

(2) W is symmetric with respect to the interchange: i-^ {v,u) 

(3) W is analytic in its variables near u — v = 0, such that W = 0(4). 

The first property is justified by the standard derivation of the coupled-mode sys- 
tem (|2.1|) with an envelope approximation The second property defines a class 
of symmetric nonlinear potentials. Although it is somewhat restrictive, symmetric 
nonlinear potentials are commonly met in physical applications of the system H2.1|l . 
The third property is related to the normal form analysis , where the nonlinear 
functions are approximated by Taylor polynomials. Since the quadratic part of the 
potential function is written in the left-hand-side of the system (|2.1ll and the cu- 
bic part violates the gauge transformation and analyticity assumptions, the Taylor 
polynomials of W start with quadric terms, denoted as 0(4). 

We find a general representation of the function W{u,u,v,v) that satisfies the 
conditions (l)-(3) and list all possible (four-parameter) quadric terms of W. 



2. Coupled-mode system 



We consider the Hamiltonian coupled-mode system in the form: 



(2.1) 
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Lemma 2.1. IfWCzC and property (1) is satisfied, such that 

(2.2) W{u,u,v,v) = W {ue'°',ue-'°',ve'°',ve-'°') , Va G M, 
then W = W{\u\^, \v\^,uv). 

Proof. By differentiating (|2.2|l in a and setting a = 0, we have the differential 
identity: 

(2.3) DW = i ( u— - u— + v—-v— \ W{u, u, v, v) = 0. 

\ ou ou ov av J 

Consider the set of quadratic variables 

which is independent for any u ^ and w 7^ in the sense that the Jacobian is non- 
zero. It is clear that 2:1.2, 3 = and DZ4 = 22:4. Therefore, DW = 2z4dz^W = 0, 
such that W = W{zi, Z2, z^). □ 

Corollary 2.2. I/WeM. and property (1) is met, then W = W{\u\'^ , uv^vu). 

Lemma 2.3. If W e M. and properies (l)-(S) are satisfied, then W — W{\u\'^ + 

Proof. By Corollary 12.21 and property (2) , we can re-order the arguments of W as 
W — W{\u\ + \v\,\u\\v\,uv + vu). By analyticity in property (3), W may depend 
only on |up and |wp rather than on |u| and \v\. □ 

Corollary 2.4. IfW G M and properties (l)-(3) are satisfied, then 

Corollary 2.5. The only quadric potential function £ M that satisfies properties 
(l)-(3) is given by 

(2.5) W = y (|u|4 + \v\^) + a2\u\^\v\^ + aad^P + \v\^){vu + vu) + ^{vu + vuf, 

where (01,02,03,04) are real-valued parameters. It follows then that 
d:^W = oilupu -I- 02-u|wp + 03 [(2|u|2 + \v\'^)v + u'^v] + 04 \v'^u + 



a^W = ai\v\^v + 02f |up + 03 [(2|w|2 + |u|2)u + y^U 



- 04 [u'^v + 



The potential function H2.5|l with 01,02 7^ and 03 = 04 = represents a 
standard coupled-mode system for a sub-harmonic resonance, e.g. in the context 
of optical gratings with constant Kerr nonlinearity pp. When 01 = 03 = 04 = 0, 
this system is integrable with inverse scattering and is referred to as the massive 
Thirring model |T2|. When oi = 02 = and 03,04 ^ 0, the coupled-mode system 
corresponds to an optical grating with varying, mean-zero Kerr nonlinearity, where 
03 is the Fourier coefficient of the resonant sub-harmonic and 04 is the Fourier 
coefficient of the non-resonant harmonic ^ (see also 11). 

We rewrite the coupled-mode system (|2.1|l as a Hamiltonian system in complex- 
valued matrix- vector notations: 

(2.6) ^ = JVi7(u), 
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where u = (u, m, v, , 





1 


i 



and (w, li, v) = Jj^ u, i;, i;)dx is the Hamiltonian functional with the density: 

z 1 
h = u, V, v) - {vu + uv) + -(uuj; - Uxu) - -(vv^ - VxV). 

The Hamiltonian H(u,u,v,v) is constant in time t > 0. Due to the gauge invari- 
ance, the coupled-mode system (|2.1(l has another constant of motion Q(u,u,v,v), 
where 

(2.7) Q = J {\u\^ + \v\^)dx. 

Conservation of Q can be checked by direct computation: 
(2.8) 



^^i\uf + \vf) + ^{\u\'-\vf)^DW^O, 



where the operator D is defined in H2.3() . Due to the translational invariance, the 
coupled- mode system (|2.1f) has yet another constant of motion P{u, u, v, v), where 



(2.9) 



P 



(uUx — UxU + VVx — VxV) dx. 



In applications, the quantities Q and P are referred to as the power and momentum 
of the coupled-mode system. 



3. Existence of gap solitons 
Stationary solutions of the coupled-mode system (|2.1|l take the form: 



(3.1) 



Ust{x,t) = uo(a; + s)e''^*+*« 



where (s,0) G arc arbitrary parameters, while the solution (uo,fo) G on 
a; S R and the domain for parameter a; G R are to be found from the nonlinear 
ODE system: 



(3.2) 



Stationary solutions are critical points of the Lyapunov functional: 

(3.3) A = M, w, w) -I- cj(5(w, u, w, w), 

such that variations of A produce the nonlinear ODE system H3.2|l . 

Lemma 3.1. Assume that there exists a decaying solution (uo,wo) of the system 
^3. i^) on a; G R. // G R satisfies properties (l)-(3), then uq = wo (module to an 
arbitrary phase). 
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Proof. It follows from the balance equation (|2.8|l for the stationary solutions 
that 

ko|'-|«oP = Co = 0, VxeR, 
where the constant Cq = is found from decaying conditions at infinity. Let us 
represent the solutions {uq,vq) in the form: 



(3.5) 



(3.6) 



such that 

iQ' - 2Q{Q' + $') = 2ljQ - 2Qe-2« + 2uodu„W{uo, uq, wq, vq) 
-iQ' - 2Q(e' - $') = 2ojQ ~ 2Qe2*® + 2vQd^,W{uQ, uq, vq, vq) 

Separating the real parts, we obtain 

Q(cos(2e) - - e' - $') = Re[uoduoW{uo,uo,vo,vo)] 
g(cos(2e) ~uj-e' + <P')^Re [vody„W{ua, uq, va, vq)] 

By Corollarv l2.4l we have = 0, such that = <i>o. □ 

Corollary 3.2. Let uq = vq. The ODE system f3.S\) reduces to the planar Hamil- 
tonian form: 

where p = 20, q — Q, and 

(3.8) h = W{p,q)-2qcosp + 2ujq, W{p, q) ^ W{uo,uq,vo,vo). 

Proof. In variables {Q, &) defined by H3.4|l with <I>(a;) = $o = 0, we rewrite the 
ODE system (|3.5ll as follows: 

Q' = 2Q sin(2e) + 2Im [uoduoW{uo, uo, vo, vo)] 
QQ' = ~ujQ + Qcos(20) - Re[uoduoW {uo,uo,vo,vo)] 

The system (|3.9|) is equivalent to the Hamiltonian system H3.7|l and H3.8|l if 

dpW{p, q) = i [uodua - uoduo] W{uo, uq, vo,vq) 
qdqW{p, q) = [uaduo + uoSoo] W{uo, uo, vq, vq) 
The latter equations follows from H2.3|l . (|2.4|l . and 13. 4() with the chain rule. □ 
Corollary 3.3. Let uq ~ vq. Then, 

(3 11) d'^ - W = d'^ - W d'^iW = d'^^W 9^ W = d"^ - W 

Remark 3.4. The family of stationary solutions (|3.1|l can be extended to the family 
of travelling solutions of the coupled-mode system H2.1|l by means of the Lorcntz 
transformation ^S]. With the boosted variables, 

Ct „ t — CX /I— C\ I 1 + 



(3.9) 



(3.10) 



X^^==, T=^==, U=^^ u, 



where c G (—1,1), the family of travelling solutions still satisfies the constraint 
\Uo\'^ = iKjP from the balance equation (|2.8|) . However, CoroUarv 12.41 fails for a 
boosted potential function W{U, U, V, V) and the representation H3.4|l results no 
longer in the relation J7o = El- It will be studied separately if the block- 
diagonalization of the linearized coupled-mode system can be extended (in a non- 
trivial matter) to the family of travelling solutions. 
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Decaying solutions of the system (13.2(1 with a homogeneous polynomial function 
W{u, u, V, v) are analyzed in Appendix A. Conditions for their existence are identi- 
fied for the quadratic potential function (|2.5|l . Decaying solutions may exist in the 
gap of continuous spectrum of the coupled- mode system (|2.1|l for cu G (—1,1). We 
introduce two auxiliary parameters: 

l-UJ 



(3.12) A*=T— ' P=VT^, 

such that < ^ < oo and < /3 < 1. When ai ~ 1, a2 — p, and — — 0, we 
obtain in Appendix A the decaying solution uo{x) in the explicit form: 



I2{1-LU) 1 



'^^ \l I + p {cosh (3x + iy/JI sinh (3x) 

When u! ^ 1 (such that p, ^ and f3 ^ 0), the decaying solution (|3.13() becomes 
small in absolute value and approaches the limit of sech-solutions sech(/3x) . When 
to —> —1 (such that /i — > oo and /3 — > 0), the decaying solution H3.13|l remains finite 
in absolute value and approaches the limit of the algebraically decaying solution: 

2 

Un = — , -. 

VTTp(l + 2ix) 

When oi = a2 = 0, 03 = 1 and 04 — s, the decaying solution 1*0(2;) exists in two 
sub-domains: w > 0, s > — 1 and cj < 0, s < 1. When w > 0, s > —1, the solution 
takes the form: 



^ 1 — a; (cosh /3x — 1^/1 sinh /3x) 



2 

where 

A+ = [(s - l)//2 _ 2sp +{s + 1)] cosh^(/3x) + 2[sp - (s - l)p^] cosh^ (/3a;) + (s - l)p^. 
When u! < 0, s < 1, the solution takes the form: 

ll — UJ (sinh /3a; — i^/Ecosh Bx) 

'"^> "° ^ V— ^ 

where 

A_ = [(s + 1) - 2sp - (s - l)p^] cosh''(/3a;) + 2[s + l- sp] cosh^(/3a;) - (s + 1). 

In both limits w — > 1 and a; — 1, the decaying solutions H3.14|) and H3.15|l approach 
the small-amplitude sech-solution sech(/3a;). In the limit — s- 0, the decaying 
solutions (|3.14(l and H3.15|l degenerate into a non-decaying bounded solution with 
K(a;)P = i. 

4. BLOCK-DIAGONALIZATION of the LINEARIZED SYSTEM 

Linearization of the coupled-mode system H2.1|l at the stationary solutions (|3.1|) 
with s = 6* = is defined as follows: 

u{x, t) = e^'^* {uq{x) -f Ui{x)e^^) 
I u{x,t)^e--^*{uo{x) + U2{x)e^*) 
^ ' ^ v{x, t) = e"^* {vq{x) + U3{x)e^*) 

u(x,t) = e~"^* {vo{x) + U4ix)e^*) 
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where vq — uq, according to Lemma l3. II Let (f , g) be a standard inner product for 
f , g G L^(]R,C'*). Expanding the Lyapunov functional into Taylor series near 
uo = (uq, uo, t;o, wo)^, we have: 

1 



(4.2) 



A = A(uo) + (U,VA|„J 



(U,i/^U) 



where U = {Ui,U2, 113,1/4)'^ and is the the linearized energy operator in the 
explicit form 



(4.3) 
where 

(4.4) 
and 

(4.5) 



D = 



Lo — idx 


-1 






uj + idx 


-1 



-1 


LJ + idx 




idx ) 



V = 



UqUo 
VQUO 

,2 

VoUo 



V 9- 



d\ 



92 
92 
92 

d\ 



92 - 
92 - 



92 - 



l¥(uo,Uo,Wo,Wo)- 



The linearization 14. 1|) of the nonlinear coupled-mode system (|2.1|l results in the 
linearized coupled-mode system in the form: 



(4.6) 



where cr is a diagonal matrix of (1,-1,1,-1). Due to the gauge and translational 
symmetries, the energy operator i?^ has a non-empty kernel which includes two 
eigenvectors: 



(4.7) 



Ui (TUo(a;), U2 = UQ(a;). 



The eigenvectors U1.2 represent derivatives of the stationary solutions H3.1|) with 
respect to parameters (6*, s). We adopt a standard assumption that the coupled- 
mode system is generic. 

Assumption 4.1. The kernel of is exactly two-dimensional with the eigenvec- 
tors u^ . 

Due to the Hamiltonian structure, the linearized operator aHi^ has at least 
four-dimensional generalized kernel with the eigenvectors (|4.7|l and two generalized 
eigenvectors (see IS"* for details). The eigenvectors of the linearized operator aH^ 
satisfy the tr-orthogonality constraints: 



(4.8) 
(4.9) 



(uo,U) 
K,f^U) 



{uqUi + U0U2 + V0U3 + V0U4) dx = 0, 
{uqUi — u'qU2 + VqUs — v'qU4) dx = 0. 



The constraints 14.8|l and (|4.9|l represent first variations of the conserved quantities 
Q and P in (j2.7|l and H2.9|l at the linearization H4.1|l . 
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It follows from the explicit form of and from Corollary that the eigenvalue 
problem H^oXJ = /iU has two reductions: 

(4.10) {i) Ui = Ui, U2 = U3, {li) Ui - -Ui, U2 - -Us. 

Our main result on the block-diagonalization of the energy operator H^^ and the 
linearized coupled- mode system ^AM is based on the reductions (|4.1Q(I . 

Theorem 4.2. Let € M satisfy properties (l)-(3). Let (uo,wo) &e a decaying 
solution of the system ^3.^) on x € M, where vq = Uq- There exists an orthogonal 
similarity transformation S, such that S^^ — , where 

/ 1 1 \ 
10 1 
10-1 

V 1 -1 y 

that simultaneously block- diagonalizes the energy operator H^, 



V2 



(4.11) S-'H^S ^y^Q ]^H, 
and the linearized operator aHi^ 

(4.12) S-'aH^S = <j( ) = iL, 



, H+ 

where H± are two-by-two Dirac operators: 



an- 



d 



pp. -1-^2 a2 _|_ 32 \ 

(4.14) V± = pa j^pa pa° j^pfl W^(uo, uq, wo, wo)- 

\ '-'ul ^ "«oSo '^"o"o ^ "novo J 

Proof. Applying the similarity transformation to the operator D{dx) in (|4.4|l . we 
have the first terms in Dirac operators H±. Applying the same transformation to 
the potential V{x) in H4.5|l and using Corollarv l3.3l we have the second term in the 
Dirac operators H± . The same transformation is applied similarly to the linearized 
operator aLl^^ with the result H4.12|l . □ 

Corollary 4.3. The linearized coupled-mode system \4.()\ is equivalent to the block- 
diagonalized eigenvalue problems 

(4.15) a^H^asH+Vi = 7V1, (jgiJ+fTa-ff- V2 = 7V2, 7 = -A', 
where Vi^2 G and is the PaulVs diagonal matrix of (1, —1). 

Corollary 4.4. Let Uo = (iio,So) G and (f,g) be a standard inner product for 
f , g G L^(M,C^). Dirac operators H± have simple kernels with the eigenvectors 

(4.16) H+u'„ = 0, i/-CT3Uo = 0, 
while the vectors V1.2 satisfy the constraints 

(4.17) (uo,Vi)=0, (u;„a3V2) = 0. 
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Remark 4.5. Block-diagonalization described in Theorem 14.21 has nothing in com- 
mon with the exphcit diagonahzation used in reduction (9.2) of 14 for the partic- 
ular potential function H2.5|l with ai ~ a2 — = and ~ 1. Moreover, the 
reduction (9.2) of ^l] does not work for oj ^ 0, while gap solitons do not exist in 
this particular model for uj — 0. 

We illustrate block-diagonalization of the eigenvalue problem (14.151) for the quadric 
potential function H2.5(l . When ai — 1, a2 — p and 03 = 04 = 0, the decaying solu- 
tion uq{x) is given by H3.13|l and the potential matrices V±(x) in the Dirac operators 
H± in H4.13|l - I|4.14(l are found in the explicit form: 

When oi = 02 = 0, 03 = 1 and 04 = s, the decaying solution ui^{x) is given by 
either 13.1411 or H3.15() and the potential matrices V±(x) take the form: 



(4.19) = 3 



(4.20) V- = 



-2|?/o|' 



^0 



Numerical computations of eigenvalues of the Dirac operators Ii± and the linearized 
operator L in l|4.11(l and (|4.12(l are developed for the explicit examples (|4.18(l and 

5. Numerical computations of eigenvalues 

Numerical discretization and truncation of the linearized coupled-mode system 
(|4.6I) leads to an eigenvalue problem for large matrices J^I- Parallel software li- 
braries were recently developed for computations of large eigenvalue problems |20j . 
We shall use Scalapack library and distribute computations of eigenvalues of the 
system ()4.6|l for different parameter values between parallel processors of the SHAR- 
Cnet cluster Idra using Message Passing Interface (211 ■ 

We implement a numerical discretization of the linearized coupled-mode system 
(|4.t)|) using the Chebyshev interpolation method [22]. The main advantage of the 
Chebyshev grid is that clustering of the grid points occurs near the end points 
of the interval and this clustering prevents the appearance of spurious complex 
eigenvalues from the discretization of the continuous spectrum. If the eigenvector 
is analytic in a strip near the interpolation interval, the corresponding Chebyshev 
spectral derivatives converge geometrically, with an asymptotic convergence factor 
determined by the size of the largest ellipse in the domain of analyticity ^| . 

The continuous spectrum for the linearized coupled-mode system (|4.6|) can be 
found from the no-potential case Vix) = 0. It consists of two pairs of symmetric 
branches on the imaginary axis X G iM. for |Im(A)| > 1 — uj and |Im(A)| > 1 + uj 
[71E1- In the potential case V{x) ^ 0, the continuous spectrum does not move, but 
the discrete spectrum appears. The discrete spectrum is represented by symmetric 
pairs or quartets of isolated non-zero eigenvalues and zero eigenvalue of algebraic 
multiplicity four for the generalized kernel of aH^^ [ZIE]- We note that symmetries 
of the Chebyshev grid preserve symmetries of the linearized coupled-mode system 
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We shall study eigenvalues of the energy operator Hui , in connection to eigenval- 
ues of the linearized operator aH^j- It is well known ^1|25] that Hermitian matrices 
have condition number one, while non- Hermitian matrices may have large condition 
number. As a result, numerical computations for eigenvalues and eigenvectors have 
better accuracy and faster convergence for self-adjoint operators ^1|22]- We will 
use the block-diagonalizations H4.11|l and H4.12|l and compute eigenvalues of 

and L. The block-diagonalized matrix can be stored in a special format which 
requires twice less memory than a full matrix and it accelerates computations of 
eigenvalues approximately in two times. 

Figure 1 displays the pattern of eigenvalues and instability bifurcations for the 
symmetric quadric potential (|2.5(l with oi = 1 and 02 = 03 = a4 = 0. The decaying 
solution ^0(2:) and the potential matrices V±{x) are given by H3.13() and (|4.18(l 
with p = 0. Parameter oj of the decaying solution ^0(2;) is defined in the interval 
— 1 < cij < 1. Six pictures of Fig. 1 shows the entire spectrum of L, and 
H- for different values of u. (The continuous movie that shows transformation 
of eigenvalues when uj decreases is available as a multi-media attachment to this 
article.) 

When oj is close to 1 (the gap soliton is close to a small-amplitude sech-soliton), 
there exists a single non-zero eigenvalue for and and a single pair of purely 
imaginary eigenvalues of L (see subplot (1) on Fig. 1). The first set of arrays on 
the subplot (1) indicates that the pair of eigenvalues of L becomes visible at the 
same value of uj as the eigenvalue of H^. This correlation between eigenvalues of 
L and can be traced throughout the entire parameter domain on the subplots 
(1H6). 

When bj decreases, the operator acquires another non-zero eigenvalue by 
means of the edge bifurcation ^3], with no changes in the number of isolated eigen- 
values of L (see subplot (2)). The first complex instability occurs near oj ~ —0.18, 
when the pair of purely imaginary eigenvalues of L collides with the continuous 
spectrum and emerge as a quartet of complex eigenvalues, with no changes in the 
number of isolated eigenvalues for and H- (see subplot (3)). 

The second complex instability occurs at w « —0.54, when the operator 
acquires a third non-zero eigenvalue and the linearized operator L acquires another 
quartet of complex eigenvalues (see subplot (4)). The second set of arrays on the 
subplots (4)-(6) indicates a correlation between these eigenvalues of L and 

When bj decreases further, the operators and acquires one more isolated 
eigenvalue, with no change in the spectrum of L (see subplot (5)). Finally, when to 
is close to — 1 (the gap soliton is close to the large-amplitude algebraic soliton) , the 
third complex instability occurs, correlated with another edge bifurcation in the 
operator (see subplot (6)). The third set of arrays on subplot (6) indicates this 
correlation. The third complex instability was missed in the previous numerical 
studies of the same system [71E|. In a narrow domain near ui ~ —1, the operator 
has two non-zero eigenvalues, the operator i7_ has five non-zero eigenvalues 
and the operator L has three quartets of complex eigenvalues. 

Figure 2 displays the pattern of eigenvalues and instability bifurcations for the 
symmetric quadric potential (|2.5(l with ai = 02 = 04 = and 03 = 1 . The decaying 
solution uq (x) and the potential matrices V± (x) are given by H3.14|) and H4.19|l with 
w > and s = 0. Eigenvalues in the other case lo < can be found from those in 
the case w > by reflections. 
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When Lo is close to 1 (the gap sohton is close to a small-amplitude sech-soliton), 
there exists one non-zero eigenvalue of and no non-zero eigenvalues of L and 
(see subplot (1)). When oj decreases, two more non-zero eigenvalues bifurcate 
in H- from the left and right branches of the continuous spectrum, with no change 
in non-zero eigenvalues of L (see subplot (2)). The first complex bifurcation occurs 
at w « 0.45, when a quartet of complex eigenvalues occurs in L, in correlation 
with two symmetric edge bifurcations of from the left and right branches of 
the continuous spectrum (see subplot (3)). The first and only set of arrays on the 
subplots (3)-(6) indicates a correlation between eigenvalues of L and H^, which 
is traced through the remaining parameter domain of lo. The inverse complex 
bifurcation occurs at uj « 0.15, when the quartet of complex eigenvalues merge at 
the edge of the continuous spectrum into a pair of purely imaginary eigenvalues 
(see subplot (5)). No new eigenvalue emerge for smaller values of lo. When lo is 
close to (the gap soliton is close to the non-decaying solution), the operator 
has two non-zero eigenvalues, the operator H- has three non-zero eigenvalues and 
the operator L has one pair of purely imaginary eigenvalues (see subplot (6)). 

We mention two other limiting cases of the symmetric quadric potential 1)2. 5|l . 
When oi = 03 = a4 = and a2 = 1, the coupled-mode system is an integrable 
model and no non-zero eigenvalues of L exist, according to the exact solution of the 
linearization problem pi llOj . When oi = 02 = 03 = and 04 — ±1, one branch of 
decaying solutions ^0(2:) exists for either sign, according to (|3.14() and H3.15|l . The 
pattern of eigenvalues and instability bifurcations repeats that of Fig. 2. 

Numerical results reported above imply that the number of isolated non-zero 
eigenvalues of the linearized operator L is bounded from above by the total number 
of non-zero isolated eigenvalues of the energy operators H-^- and H^. Furthermore, 
there exists a correlation between edge bifurcations in the operator L and those in 
the Dirac operators and These analytical questions will be addressed in 
the future work. 



Appendix A. Conditions for existence of gap solitons in the 

HOMOGENEOUS POTENTIAL FUNCTION 

We shall consider the homogeneous potential function g R of the monomial 
order 2n that satisfies properties (l)-(3). The general representation of W{u, u, v, v) 
is given by 



(A.l) W^ = ^^afe,,(?/^«^ + S^z;^)|M| 



2n-2k-2s\y\2k 



s=0 fc=0 



where ak,s are real- valued coefficients which are subject to the symmetry conditions: 
o-ki,s — cik2.s ii ki+k2 = s for 3 = 0,1,.. .,71— 1. Let's introduce new parameters 



n — s 



As^^ak,s, s = 0,l,...,n. 



fe=0 



Using the variables (Q,6) defined in (|3.4() with $(a;) = $0 = 0, we rewrite the 
ODE system 1)3. 7|l in the explicit form: 



(A.2) 



Q' = 2Q sin(2e) - 2Q" X^Lo sin(2se) 
e' ^-Lo + cos(2e) - nQ"-i X;"=o cos(2se) 
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There exists a first integral of the systera (|A.2|I : 

n 

-uiQ + cos(2e)Q ~Q"^As cos(2se) = Co, 

s=0 

where Co = from the zero boundary conditions Q{x) — s- as |a;| — > oo. As a 
result, the second-order system (|A.2ll is reduced to the first-order ODE 

(A.3) e'(j:) = (n- l)(w-cos(2e)), 

while the function Q{x) > can be found from Q{x) as follows: 

We consider the quadric potential function W given by H2.5() . Using ljA.3|) for 
the case n = 2 we obtain: 

(A.5) e'(a;) = to - cos(2e), 

and the correspondence: 

ai+a2+ aj 4 "4 

^0 = ^ , Ai = 2a3, A2 = —. 

We rewrite the representation (|A.4p for Q{x) as follows: 
where 

t = cos(2e), = a4t2-|-2a3t+^^i^i-^, 

such that i G [—1,1]. Let's consider two cases: 

^ 1 i <w; 0(0 <0 =>Q- 

We can solve the first-order ODE (|A.5p using the substitution z — tan(8), such 
that 

t = ^ = . 

1 + z2 1 + i 

After integration with the symmetry constraint 0(0) — 0, we obtain the solution 

(^- Vm) 



(A.8) 
where 



/5 = v/r^, ^, 



l+UJ 

and — 1 < < 1. Two separate cases are considered: 

(A.9) |z|<Vm - = t= ^"^^;(^")-^""^;(^"\ 

^^cosh(/3a:;) cosh^ (/Jx) -f ^ sinh^ (/3a;) 

where t > lu, and 

CA inA 11-.^ - msh(£x) ^ _ smh^jPx) ~ ficosh^{(3x) 
(A.IO) |z|>Vm ^"sinh2(/3a;)+/zcosh^(/3a;)^ 
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where t < u. Let's introduce new parameters 

ai + 02 



A = — 2a3 + 04 
B = — 2a4 + a 
C = 2a3 + tti 



2 

B = -2a4 + fli + a2, 



2 

It is clear that A — 0(— 1) and C = '/'(I)- If i > w and (/)(t) > 0, it follows from 
(TO)! and (fO)l that 

(A.ll) Q+{x) = 



{An^ + Bfi + C) cosh''(/3x) - {Bfi + 2An^) cosh2(/3a;) + A^^ ' 
Et<uj and < 0, it follows from (fO)l and (|A.10|I that 

M 12) O-fx) = i^-mt^+l)cosh\(3x)-l) 

^ ' ' ^ ' {Afj.'^ + Bn + C)cosh'^{(3x) - {Bfi + 2C)cosh^if3x) + C' 

The asymptotic behavior of the function Q{x) at infinity depends on the location 
of the zeros of the function V'(m) = ^/^^ + B11 + C. The function V"!/^) is related to 
the function (t>{t), e.g. if = then 4'{u!) — 0. 

A.l. Case A < 0, C > 0. In this case the quadratic polynomial (p{t) has exactly 
one root 4'{ti) — such that ti G (—1,1). We have two branches of decaying 
solutions with the positive amplitude Q{x). One branch occurs for ti < lu < 1 with 
Q{x) = Q'^{x) and the other one occurs for —l<uj<ti with Q{x) — Q^{x). At 
the point oj = ti, the solution is bounded and decaying. 

A. 2. Case yl > 0, C > 0. In this case the quadratic polynomial (f){t) has no roots 
or has exactly two roots on (—1, 1). If does not have any roots on (—1, 1), we 
have a decaying solution with the positive amplitude Q{x) for any — 1 < w < 1 
with Q{x) = Q^{x). If has two roots (/'(ti) — and 0(i2) = such that 
^1,^2 G 1) then we have a decaying solution with Q{x) ~ Q^{x) only on the 
interval max(ti,t2) < uj < 1. At the point lo — max(ti,i2), the solution becomes 
bounded but non-decaying if ti ^ t2 and unbounded if ti — t2. 

A. 3. Case A < 0, C < 0. In this case the quadratic polynomial (f>{t) has no roots 
or has exactly two roots on (—1, 1). If does not have any roots on (—1, 1), we 
have a decaying solution with the positive amplitude Q{x) for any — 1 < w < 1 
with Q{x) = Q^{x). If (t>{t) has two roots (t>{ti) — and 0(^2) — such that 
^17^2 G (^1? 1) then we have a decaying solution with Q{x) = Q^{x) only on the 
interval —1 < a; < min(ti,i2)- At the point lo — min(ii,<2)7 the solution becomes 
bounded but non-decaying if ti ^ t2 and unbounded if ti = t2. 

A. 4. Case yl > 0, C < 0. In this case no decaying solutions with positive ampli- 
tude Q{x) exist. 

A. 5. Special cases. Two special cases occur when 0(1) = or 0(— 1) = 0. If 
<p{l) = then Q^{x) has a singularity at x = for any — 1 < < 1. If (/'(—I) = 
then Q~{x) has a singularity at a; = for any — 1 < w < 1. 
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Figure 1. Eigenvalues and instability bifurcations for the sym- 
metric quadric potential l|2.5|) with ai = 1 and 02 = 03 = 04 = 0. 
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Figure 2. Eigenvalues and instability bifurcations for the sym- 
metric quadric potential 12. 5|) with — 1 and ai ~ a2 — = 0. 
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